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Algebraic multigrid (AMG) is a numerical method used to solve particular 
algebraic systems, and interest in it has risen because of its multigrid-like efficiency. 
Variations in methodology during the interpolation phase result in differing conver- 
gence rates. We have found that regular interpolation weight definitions are inade- 
quate for solving certain discretized systems so an iterative approach to determine 
the weights will prove useful. This iterative weight definition must balance the re- 
quirement of keeping the interpolatory set of points “small” in order to reduce solver 
complexity while maintaining accurate interpolation to correctly represent the coarse- 
grid function on the fine grid. Furthermore, the weight definition process must be 
efficient enough to reduce setup phase costs. 

We present systems involving matrices where this iterative method signifi- 
cantly outperforms regular AMG weight definitions. Experimental results show that 
the iterative weight definition does not improve the convergence rate over standard 
AMG when applied to M-matrices; however, the improvement becomes significant 
when solving certain types of complicated, non-standard algebraic equations gener- 
ated by irregular operators. 
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I. 



INTRODUCTION 



Precursors to multigrid (MG) methods were developed in the mid 1960’s by 
Fedorenko (1964), and Bachvalov (1966), and remained fairly unknown until Achi 
Brandt [Ref. 1] fully developed multigrid in his 1973 article. In that paper, he 
demonstrated for the first time the utility of true multigrid methods as a fast numeri- 
cal solvers for boundary value problems. An explosion of papers ensued and the 1980’s 
produced a plethora of articles as well as efficient and reliable computer algorithms to 
further demonstrate the practicality and usefulness of MG. Confidence in multigrid 
grew when the true merits of this solver came to bear fruit in faster convergence rates 
over conventional methods. Proofs of these faster convergence rates coupled with a 
sufficient number of satisfactory numerical results gave way to a general acceptance 
of the method even to the most sceptic [Ref. 2]. 

Originally developed to solve simple boundary value problems, multigrid meth- 
ods gained wide recognition for their speed and efficiency in solving general linear 
partial differential equations (PDEs) of the elliptic form, e.g., (in two dimensions) 

-V 2 u(x,y) = f(x,y), 



or explicitly as 

^XX 

If this equation is now discretized, for example by finite differences, on some rectan- 
gular domain, ft C as in Figure 1, and if we denote the boundary of the domain 
as 5ft, then multigrid will solve problems involving Laplace’s equation 

— V 2 u = 0 in ft 
u = 0 on 5ft 



(i.i) 



i 



8 




Q I I I I I I I 

012345678 

x 



Figure 1. Rectangular domain for discretization of a 2-D problem. 



Poisson’s equation 



Helmholtz’s equation 



— V 2 t£ = S in Cl 
u = 0 on dfl 



— V 2 u + k 2 u = S in O 
u = 0 on dfl 
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or even the anisotropic Poisson equation 



d 2 u 

£ d. x* 



Pu 

dy 2 



= S in Q 



u = 0 on dQ 



( 1 . 2 ) 



(1.3) 



(1.4) 



in any number of dimensions almost effortlessly (comparatively speaking, of course). 

These types of equations frequently arise in physical applications such as 
steady state temperature problems, fluid flow, orbital mechanics (gravitational fields), 
steady state electrical field problems, and many others. When a matrix equation arises 
from the discretization of (1.1) through (1.4), multigrid methods are quite successful. 

MG methods are fast iterative solvers using a hierarchy of levels, or grids, and 
may be used with many types of discretization techniques such as finite difference 
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or finite element methods. In solving problems of the elliptic nature using these 
discretization techniques, multigrid has proven to be the fastest numerical solution 
technique in the field. Unlike other numerical solvers, multigrid is general enough so 
that it can effectively use fairly arbitrary regions and boundary conditions and more 
importantly, does not depend on the separability of the differential equations. 

In fact, from presentations at the Copper Mountain Conference on Multigrid 
Methods and from various papers which can be found on MG Net (“http://casper.cs. 
yale.edu / mgnet / www/mgnet-papers.html# Y”), we know that multigrid can also be 
directly applied to more complicated, non-symmetric and nonlinear systems of equa- 
tions, like the Lame-System of elasticity or the Stokes (or Navier-Stokes) equations. 
Multigrid has been applied successfully to electrostatic and magnetostatic problems 
[Ref. 3], to statistical physics problems, to integral equations and to image recon- 
struction algorithms [Ref. 4, 5]. It can also be applied to problems in control theory, 
partical physics and permeable magnetic materials [Ref. 2]. 

The main concept of multigrid methods is “to complement the local exchange 
of information in point-wise iterative methods (on each level) by a global one utilizing 
several related systems, called course levels , with a smaller number of variables [Ref. 
6].” The key to multigrid’s performance then is to apply a relaxation technique 
as many times as needed to dampen the oscillatory modes of the error (since we 
know relaxation works best on highly oscillatory modes) and then apply a coarse-grid 
correction scheme to eliminate the smooth components of the error. This combination 
of relaxation and coarse-grid correction is the essence of why MG works so well. 

Using multigrid methods, we can solve a large, difficult problem by reduc- 
ing it iteratively into successive smaller, easier ones. At the coarsest level then the 
system can be solved directly using a known, efficient direct method (such as the 
LU-decomposition, Gaussian elimination, etc.). For systems such as (1.1) through 
(1.4), very efficient methods are available. From A Multigrid Tutorial [Ref. 7], we 
know that when applied to an N x N grid, multigrid methods are nearly optimal 
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since they only require 0(N 2 log N) arithmetic operations and hence approach the 
minimum operation count of 0(N 2 ) operations. In fact, it is shown in [Ref. 7] that 
the full multigrid V-cycle (FMV) method which we will discuss in the next chapter, 
requires only 0(N d ) to the level of discretization on the “model” problem (Laplace’s 
equation (1.1)) of dimension d. This operation count is optimal. 
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II. THE FUNDAMENTALS OF 
MULTIGRID 

There are two fundamental components of multigrid; one is the idea of coarse- 
grid correction (CG) and the other is that of nested iteration (NI). Presented here is a 
basic outline of both methods. To begin, we note that the principle concept of CG is 
that it is a correction strategy that transfers the components of the problem between 
the fine and coarse grids. Complementing coarse-grid correction is nested iteration 
which is based on the following concept: use the information on the coarse grids to 
provide an informed guess for the initial guess on the finer level. 

A. THE PROBLEM STATEMENT 

In order to understand the basics of multigrid and to motivate its usage, let 
us apply the method to a simple one-dimensional problem. We begin by discretizing 
a problem of the form 

cPu(x) . . ... 

fa? +au W = fW’ 

with Dirichlet boundary conditions (u(a:) = 
in Figure 2, into N + 1 points: 

x j = jh, where j = 0, . . . , N. 

We make h = 1/iV of constant width throughout the interval. This creates a grid 



x=0 x=l 




Figure 2. Discretized domain of problem II.l. 



of iV — 1 interior points. At each of these interior points, the equation in (II.l) is 



0 < 2 < 1, <7>0 (II.l) 

0 on dLl) by partitioning the domain, as 
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replaced by a second order discretization, such as a finite difference approximation 



-Vj - 1 + 2 Vj - Vj + 1 

h 2 



+ <roj = fj, 1 < j < N - 1 

Vo = Vn = 0 



(II.2) 



where we define fj = f(xj). The approximation (II. 2) of (II. 1) leads to a system 
of linear equations. In compact form, it is written as Av = /, or in explicit matrix 
notation, it is written 



h 2 



2 + ah 2 —1 




Vl 




h 


— 1 2 + ah 2 —1 




V2 




h 


-1 2 + ah 2 -1 










— 1 2 + ah 2 




vn- i 




Sn-\ 



where the matrix A is tridiagonal, symmetric, positive definite and has dimension 
N — 1 x iV — 1. u is a column vector containing the approximate discretized solution 
of u(x) and /,• = /(x,) is the discretized right-hand side (RHS). 



B. NOTATION 

Before we continue with the outline it is instructive to define a few terms. Let 
$l h denote the original grid or the domain on which the original problem is defined. 
We call ft h the finest grid. The symbol, fl 2h , is defined to mean the next coarser grid 
(usually of size N/2 — 1 but certainly not restricted to that). Similarly, v h will be 
an approximation to the exact solution, u h , on £l h . This means that v 2h represents 
the approximate solution to u 2h on the next coarser level Cl 2h . This notation will be 
consistent for all such grids, £l kh . The algebraic error, e kh , on f l kh ( k = 2 n , n = 
0, 1,2, 3, . . .) is given then by e kh = u kh — v kh where k = 1 represents the finest level. 

When we translate the vectors only between two consecutive grids, we will use 
the notation Q, h to represent the finer of the two grids and $l H to denote the coarser 
one. This makes it easier to follow the grid transfer sequence vice getting caught up 
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in exponential notation of the varying grids. Now that we have dispensed with some 
of the notational formalities, let us move on to more of the important issues. 

C. ITERATIVE METHODS 

We introduce iterative methods here to illustrate the need for coarse-grid cor- 
rection. First of all, iterative methods are needed and in fact, are required since the 
use of direct methods can be impractical if the matrix A is large and sparse. The 
sought-after factors of A using these direct methods can, and tend to, be very dense. 
Even when the matrix is banded (but still sparse), algorithms that are ideal for fac- 
toring such problems may be difficult to implement [Ref. 8]. Iterative methods by 
contrast “generate a sequence of approximate solutions and essentially involve the 
matrix A only in the context of matrix- vector multiplication [Ref. 8].” Iterative 
methods are attractive and easy to use because of their simplicity in implementation 
but may be prohibitively slow when the error is smooth. 

What does it mean when we say, “when the error is smooth?” To answer this 
and to understand the need for coarse-grid correction, we must first show the power 
of iterative methods and then focus on their limitations. To avoid any confusion, we 
note here that the phrase “iterative method” is isomorphic to “relaxation scheme” 
or just simply, “relaxation.” We demonstrate relaxation with the weighted Jacobi 
iteration since it is simple in nature and offers the same benefits as other iterative 
methods; in addition, it shows the common defects of relaxation in general. Other 
iterative methods include (regular) Jacobi, Gauss-Seidel, successive over-relaxation 
(SOR) and Chebyshev semi-iterative; of course, there are many others. 
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1. The Weighted Jacobi Iteration on the “Model” 
Problem 

To begin, let us solve the following one-dimensional “model” problem which 
is problem (II. 2) with a and / set to zero. The problem then becomes 



— Uj_! + 2uj — Uj+ 1 =0, l<j<N-l 

Uq = U/V = 0. 



(II.3) 



We solve (II. 3) by using the weighted Jacobi iteration. This iteration is computed 
using 



(new) 

V 



= (1 - a?)uj old ) + UJVj, 



1 < j < N - 1 



where 



» 1 ( (old) . (old) . ,2 r\ 

= o H-l + Vj+l + k fj). 



1 < j < N - 1 



and u € 'R- is a weighting factor chosen such that 0 < u < 1. In our use of the 
weighted Jacobi iteration, we set u> = 2/3. It turns out to be the optimal weighting 
factor. 



For the moment, let us assume that we have found an approximation, v ^ 
0, to the discretized solution, u = 0, on our model problem. Since the system is 
homogeneous and linear, we know the error, e, exactly, namely — v since e = u — v = 
—v. In reality, the error will consist of many different frequencies but for simplicity, 
let us assume that e consists of only three frequency modes: one high, one low and a 
third mode in the middle. Our assumption using the three modes here is merely to 
amplify the benefits of relaxation while simultaneously exploiting the inherent defects. 



2. Frequency Modes of the Error 

Consider an initial approximation to the solution consisting of the Fourier 



modes 

( jkir 

If 

where 0 < j < N and 1 < k < N — 1. The number k represents the wave number 
of the Fourier mode. The Arth mode consists of k/2 full sine waves on the interval. 
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As we will see shortly, the wave number k (frequency mode) will play a major role 
in the overall effect of relaxation. Figure 3 illustrates the modes Vj = sin with 

wave-numbers k = 1 , 6, 12 on a grid with N = 64. Wave number k = 1 represents 
the lowest frequency while k = 63 is the highest one. We now iterate (or “relax”) on 




Figure 3. Fourier modes Vj = sin with wave-numbers k = 1 , 6, 12 on a grid 

with N — 64. The kth mode consists of k/ 2 full sine waves on the interval. i>i is 
represented by the solid line, v&, the dashed and V 12 , the dot-dashed. 

problem (II. 3) by applying the weighted Jacobi iteration with u> = 2/3 on a grid of 
size N = 64. 

But wait! Maybe we are being a bit too hasty. Let’s step back from solving the 
problem for a moment. Before we see what happens with our mixed-frequency initial 
guess, it is instructive to first watch what happens to the individual components in 
our initial guess. Figure 4 shows the application of the weighted Jacobi iteration 
metho on problem (II.3) with Vi, i> 3 , i>6 , and v \2 as separate initial guesses. 
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Figure 4. The weighted Jacobi iteration with u = 2/3 applied to the 1-D “model” 
problem with N = 64 and with initial guesses v\, v 3 , vq, v 12 . ||e||oo is plotted against 
the iteration number for 50 iterations. 



Notice that the error quickly decreases to zero within a few iterations for the 
high- frequency mode (U 12 ) and although not quite as rapidly, the error also decreases 
for the middle-frequency mode (t> 6 ) but it doesn’t quite go to zero. In extreme con- 
trast, the low-frequency modes (tq and u 3 ) appear to be relatively untouched by the 
relaxation scheme. In fact, a prohibitively large number of iterations will still not 
effectively reduce the error. It is the wonderful property of iterative methods which 
allow for the quick elimination of high-frequency modes. On the other hand, they are 
quite defective in decreasing the error when it consists solely of low-frequency modes. 

We now return to problem (II. 3) with our mixed-frequency initial guess. It is 
of the form 



V] 3 



sin [ ^ | + sin 



6j 7T 

N 



+ sin 



12j7T 

N 
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Figure 5 shows how the error of a typical approximation vector decreases using re- 
laxation alone. It is not uncommon that the error decreases rapidly, usually within 
a few iterations. This is due to the quick elimination of the high-frequency modes. 
But after this initial rapid decrease, the error is reduced much more slowly. The 
main culpritsare the low-frequency modes. “The important observation is that the 
standard iterations converge very quickly as long as the error has high-frequency com- 
ponents. However, the slower elimination of the low-frequency components degrade 
the performance of the methods [Ref. 7].” 




Figure 5. The weighted Jacobi iteration with u> = 2/3 applied to problem (II. 3) with 
N = 64 and with the initial guess Vj = | |sin (p) + sin + sin • ||e||oo 

plotted against the iteration number for 50 iterations. 

What we have seen is that iterative methods work very well for the first several 
iterations but inevitably the convergence rate slows and the methods seem to stall. 
This phenomenon is due to the fact that the rapid decrease in error in the early 
stages of the method is from the swift, efficient elimination of the oscillatory modes 
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(high-frequency components). However, once the high-frequency components have 
been removed from the error, the iteration is less effective in reducing the remaining 
smooth modes (low-frequency components). 

This property of quickly removing the oscillatory modes but leaving the smooth 
modes of the error is called the smoothing property. The smoothing property of 
iterative methods (relaxation) is a serious limitation. How do we overcome this grave 
obstacle? It is in the application of coarse-grid correction to the remaining “smooth 
error” that we remedy this problem. 

D. THE METHOD OF COARSE-GRID CORRECTION 

In light of the limitations of iterative methods, we seek a post-relaxation 
scheme. The purpose of the coarse-grid correction scheme (CG) is to eliminate the 
low-frequency modes of the error once relaxation has become ineffective. The means 
by which we do this is intergrid transfers. Intergrid transfers move vectors (and the 
matrix) from the fine grid to the coarse grid and vice-versa by means of restriction 
and interpolation operators. In “moving” to a coarser grid, we observe that the 
smooth modes of the error become higher-frequency modes on this new grid. On the 
coarse-grid now, the error has oscillatory components again; but we can effectively 
remove those modes using relaxation. This two step methodology is the basis for CG. 
To understand what it means to “move” a vector between grids, we introduce a few 
required operators. 

1. The Restriction Operator 

Starting on the fine grid, Q h , we move the problem A h v h = f h to the next 
coarser grid, , by applying an intergrid transfer to A h , v h and f h . This operation 
gives us the coarse-grid equivalent to that problem but on Cl H : A H v H = f H . In doing 
this though, we must be careful to ensure that the coarse-grid problem “be consistent 
with the differential problem in the same way as the fine-grid problem [Ref. 2].” The 
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action of moving to the coarse-grid problem is performed by the restriction operator 



I? 



n 



N - 1 



If? is a linear operator from V, N ~ l to % 2 [-1 and has rank N/2—l (see [Ref. 2, 7, 9] for 
more on intergrid transfers). The most obvious of these types of restriction operators 
is injection. Injection is a linear operator which produces the coarse-grid solution 
approximation 

v H = I”v h 

where 

v? = for 1 < j < y - 1. 

In words, injection gets the value of the coarse-grid vector directly from the associated 
fine grid point. Figure 6 shows how the injection operator acts on a discretized sine 
wave when moving the vector from fl h to f \ H . 

Another type of restriction operator is full weighting. It is also a linear operator 
from to TZ~~ X with rank N/2 — 1. Like injection, it produces the coarse-grid 

solution approximation 

v H = I ”v h 

but in this case the transfer to the coarse grid is a bit different: 

v f = \ { v 2 j- 1 + 2 v{,- + v% J+1 ) for 1 < j < y - 1. 

In the case of the full weighting restriction operator, “the values of the coarse grid 
vectors are a weighted average of values at neighboring fine grid points [Ref. 7].” 
Figure 7 shows how the full weighting restriction operator acts upon a vector when 
moving the vector from the fine grid to the coarse grid. 
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Figure 6. Restriction by injection on a discretized sine wave from Cl h to Ct H . The 
circles (o) on Q h directly become the coarse-grid points on Cl H . 

2. The Interpolation Operator 

Once on the coarse grid, Tl H , and after we complete our calculations there, we 
need to be able to return to f l h . We introduce the tool that will assist us in this task: 
the interpolation operator 

l h H : — n N -\ 

The interpolation operator is a linear operator from to Tl N ~ l and has full rank. 

Ifj moves the coarse-grid vectors up one level by the calculation 

= I h H v H 



where 

Vy = vf for 0 < j < y - 1. 
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Figure 7. Restriction by full weighting on a arbitrary vector from Q h to £l H . The 
circles (o) on Vl H are the coarse-grid points and are weighted averages of neighboring 
fine grid point values. 



The interpolation operator produces a vector on the fine grid that is a weighted 
average of its adjacent points on the coarse grid. Figure 8 graphically depicts the 
action of ifj on a discretized cosine wave. 

Now that we have our intergrid transfer operators defined, it makes sense to 
talk about moving the matrix between levels. On f l h , A h = A. But what is A H 
and how do we get it? As a definition, we will call A H the f l H version of A h or, the 
coarse-grid operator. We get A H from A h by the calculation 

A h = l%A h I h H , (II. 4) 



and we recover A h from A H by the calculation 




H tH 
1 h • 
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Figure 8. Interpolation of a cosine wave from Cl H to Q, h . The stars (*) on £l h are the 
interpolated values from their adjacent coarse-grid points on Cl H . 

This is how we move the operator A on all levels. The coarse-grid operator (II. 4) is 
called the Galerkin condition. A second special and important relationship between 
the full weighting restriction operator and the linear interpolation operator is 

/^ = c(/f) r where c € K\ (II.5) 

that is, they are transposes of each other up to a (real) constant which turns out to 
be essential. (II.4) along with (II. 5) are called the variational properties. 

We note here that there is a disadvantage in using a linear interpolation op- 
erator. In fact, Wesseling explains: 

Because of the arbitrariness in choosing the direction of the diagonals of the 
interpolation triangles, there may be a loss of symmetry, that is, if the ex- 
act solution of a problem has a certain symmetry, it may happen that the 
numerical solution does not reproduce this symmetry exactly, but only with 
truncation error accuracy. [Ref. 2] 
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However, linear interpolation is much cheaper (and easier) to use than, for example, 
bilinear or quadratic interpolation. So in our work, we only use linear interpolation. 

3. The Algorithm of CG 

To begin the coarse-grid correction algorithm, we perform a relaxation scheme 
(e.g., Gauss-Seidel or Jacobi) on Cl h on the original equation A h u h = f h with an 
arbitrary initial guess v h . This has the effect of producing a much better guess 
since it removes the oscillatory components of the error right from the start. We 
now compute the residual vector, r H , on the coarser grid r H = Ik(f h — A h v h ) and 
then solve A H e H = r H for the algebraic error, e H on fl H . Now that the error e H is 
computed, we return to fl h by using the linear interpolation operator, /^, defined 
above. Since the error was smooth on Q H , I ^ transfers that error very accurately 
back to f l h . We now correct the fine grid approximation v h * — v h + e H by adding 
the interpolated error e h (= I}jC H ) back into the initial guess. Since u h = v h + e h , 
we have improved our fine grid approximation to the solution and in theory, we have 
found the solution. To improve it further, we perform relaxation on the original 
equation A h u h = f h with the updated guess v h to obtain a final approximation to the 
exact discretized solution u h . The coarse-grid correction scheme (in compact form) 
is printed below and is reproduced from Briggs’, A Multigrid Tutorial [Ref. 7 ]. 

v h CG(v h , f h ) 

Relax v\ times on A h u h — f h on f l h with initial guess v h . 

Compute r H = I?(f h - A h v h ). 

Solve A H e H = r H on fl H . 

Correct the fine grid approximation: v h < — v h + . 

Relax 1/2 times on A h u h = f h on f l h with (updated) initial guess v h . 

Here, v i and v 2 are small positive integers, usually between one and three. 
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E. THE MULTIGRID V-CYCLE SCHEME 

The multigrid V-cycle scheme is a recursive algorithm that iteratively imbeds 
the two level CG algorithm within itself. The physical property of this scheme (the 
pattern it creates) is how it gets its name. Figure 9 shows a graphic illustration of 
the V-cycle pattern. 




Figure 9. The multigrid V-cycle. 



We start one CG algorithm on the finest level and a new one on each successive 
level all the way down to the coarsest level where we the problem is solved directly. 
After we solve the problem on the coarsest level, we now complete each CG algorithm 
at each successive level all the way back to the finest level. To clarify what the V-cycle 
is really doing, we present the algorithm in a more structured format. In the code, 
it is assumed again that there are Q > 1 grids with the coarsest grid spacing given 
by Lh, where L = 2^ -1 . This procedure is also a reproduction from Briggs’ tutorial 
[Ref. 7], 
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V(v\f h ) 



Relax on A h u h = f h i/\ times with initial guess v h . 

Compute f 2h = Il h r h . 

Relax on A 2h u 2k = f 2h v x times with initial guess v 2h = 0. 

Compute f 4h = I 4 ^r 2h . 

Relax on A 4h u 4h = f 4h U\ times with initial guess v 4h = 0. 

Compute f sh = I*£r 4h . 

Solve A Lh u Lh = f Lh . 

Correct v 4h < — v 4h + I^v &h . 

Relax on A 4h u 4h — f 4h V 2 times with initial guess v 4h . 

Correct v 2h * — v 2h + 

Relax on A 2h u 2h = f 2h v <2 times with initial guess v 2h . 

Correct v h < — v h + 

Relax on A h u h = f h times with initial guess v h . 

For all its simplicity, this algorithm is the basic structure of the multigrid 
methodology. It takes individual techniques and well known concepts (relaxation and 
intergrid transfers) inherent with individual defects and integrates them in a cohesive 
way, capitalizing on the benefits of each to produce an algorithm that is extremely 
efficient and quite powerful. 

Given the recursive nature of the V-cycle algorithm above, we now present it 
in a more compact form. 
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V k +— V(v h J h ) 

1. Relax v\ times on A h u h = f h with a given initial guess v h . 

2. If Vt h = the coarsest grid, then goto 4. 

Else f 2h *- I 2 h h (f h - A h v h ) 

v 2h *- 0; 

v 2h V 2h (v 2fl , f 2fl ). 

3. Correct v h *— v h + I% h v 2h . 

4. Relax v 2 times on A h u h = f h with initial guess v h . 

The V-cycle is just one of many multigrid cycling schemes. To see its effec- 
tiveness, we show in Figure 10 how the V-cycle algorithm tailors itself to the problem 
vice the size of the matrix involved. The number of V-cycles required to solve the 
problem to truncation error barely grows even though we considerably increase the 
matrix size. 




Figure 10. Size of the problem verses number of V-cycles performed. 
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F. THE STRATEGY OF NESTED ITERATION 

When no information about the solution to the problem is known before-hand, 
it is computationally wasteful to begin solving the problem on the finest grid with just 
any initial random guess, v h . If a poor choice is made, the speed of convergence may 
be strongly affected. The algorithm may altogether fail to converge if constrained 
by the number of iterations especially when the problem is nonlinear. Computing 
the residuals on the “smaller” coarser grids is so much cheaper and it makes more 
sense than to remain only on Q h . Therefore, it is better to use the information of 
the coarse grids to provide an informed guess v h on the finest grid. In addition to 
providing a better initial guess on fl h , nested iteration gives us a better choice for v kh 
at each level k [Ref. 2]. This idea of using information inherent in the problem from 
the coarser grids to develop better initial guesses is the basis of nested iteration. In 
compact form, the NI algorithm is 

v h 4 — NI{v h J h ) 

For k = 2,4,8, 16 ,..., L 

Compute f kh = I k ( t-x)hf k ~ l)h - 

end 

Relax v\ times on Av = f on the coarsest level f l Lh and solve 

A Lh v Lh _ J-Lh for v Lh' 

Interpolate v Lh to obtain an initial guess for the Q,^ L ~^ h problem. 

For k = i — 1, ... ,8,4,2, 1 
Compute 

Relax v 2 times on Av = f on ti kh to obtain an initial guess for the Q( k ~ 1 ) h 
problem. 

end. 

In this algorithm, V\ and v 2 are again small positive integers and there are Q > 1 
grids with the coarsest grid spacing given by LA, where L = 2^ -1 . 
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G. THE FULL MULTIGRID V-CYCLE 



In this section, we present a powerful algorithm that combines nested iteration 
with the V-cycle. We also talk about how it works and why it is necessary. The 
algorithm to which we are referring is called the full multigrid V-cycle (FMV) and 
it is the method which is computationally of optimal order. We already saw the 
algorithms for NI and the V-cycle and reasoned for their need. Now we join them 
into a single, efficient method. 

In the FMV algorithm, the first approximation to the solution is obtained 
by interpolation of the solution from the next coarser grid, which has already been 
computed by another FMV algorithm. This is where NI plays its role. Once we 
achieve a first approximation, we then apply the multigrid V-cycle. “This combination 
should suffice to give a fine-grid solution to the level of truncation error [Ref. 10].” 

The FMV algorithm in compact form is (reproduced from [Ref. 7]): 

v h < — FMV h (v\ f h ) 

1. If 17^ = the coarsest grid, then goto step 3. 
else 

ph j2h(p _ A h v h ) 

v 2h «- 0; 

v 2h «- FMV 2h {v 2h , f 2h ). 

2. Correct v h «- v h + I$ h v 2h . 

3. v h <— V h (v h , f h ) vq times. 

Figure 11 shows how the FMV scheme operates over five grids. Notice that 
each V-cycle is proceeded by a smaller V-cycle which is designed to provide the best 
initial guess for the following cycle. Shortly, we will see that the extra work required in 
the preliminary V-cycle is not only of minimal cost, but it will generally pay for itself 
over the course in solving the entire problem [Ref. 7]. We also note here that FMG 
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cycles are less sensitive than the multigrid cycles but with one cycle per refinement, 
we can still guarantee a solution to the level of truncation error [Ref. 10]. 




n h 



n 2h 



n 4h 



n 8h 



n l6h 



fl 32h 

Figure 11. The full multigrid V-cycle. The (*) signify work performed before CG 
whereas (o) is the post coarse-grid correction work. The (+) on the first branch 
represent the restriction of f; no relaxation is performed here. 



Although the FMV algorithm costs more per cycle than the regular V-cycle, 
it is also more effective. “The key observation in the FMV argument is that before 
the Cl h problem is even touched, the f l 2k problem has already been solved to the level 
of truncation [Ref. 7].” The reason is that we now have a better initial guess for the 
next finer grid from our use of nested iteration. In fact, we know that because of this 
“extra” work during the preliminary cycling through the coarser grids, we only require 
0(1) V-cycles by the time we finally reach the finest grid. So the computational cost 
of the FMV method applied to a d-dimensional problem, as mentioned before, is 
0(N d ), which is optimal vice the cost of the V-cycle method alone which is of order 
0(N d log N). 
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H. IS MULTIGRID ENOUGH? 

The heart of multigrid intertwines relaxation methods such as Jacobi or Gauss- 
Seidel with intergrid operators like coarse-grid correction. In consideration of a gen- 
eral all-purpose elliptic equation solver, the two techniques enhance each other’s ef- 
fectiveness such that the overall algorithm - multigrid - is more robust and powerful 
than the sum of its parts. To be more specific, it is a known fact that relaxation 
effectively smoothes highly oscillatory modes of the error but it is ineffective on the 
smooth modes. This knowledge, coupled with the fact that any vector that lies in 
the range of the interpolation operator also lies in the null space of the coarse-grid 
correction operator, provides us a process that ensures that both the smooth and 
oscillatory modes of the error are suitably damped. 

Multigrid methods have their limitations, though. Suppose we wish to solve a 
non-elliptically structured problem or even a regular elliptic problem but defined on 
an irregular grid. How, for example, can we apply multigrid methods to these types 
of systems. Will it even work at all? Can we tailor the irregularities in the meshes? 
How do we apply MG to problems where no geometric structure (or grid) exists? 
Some specialized multigrid codes have been written for specific, ill-structured prob- 
lems to satisfy the irregular needs of the individual problem but that defeats the issue 
of providing a solver that is more nearly a “black-box” type of solver. Every special- 
ization restricts the scope of the discretization method with which the multigrid code 
can be used [Ref. 11]. A significant effort may be necessary to prevent information 
loss in the discretization and still, it may be prone to errors. It is very difficult “to 
construct an efficient preconditioning method (which modifies the condition number 
of the matrix A) for arbitrary ‘purely algebraic’ problems [Ref. 9].” 

What we really want to do is to solve purely algebraic problems with the same 
speed, efficiency and reliability that multigrid provides to geometrically structured 
problems. We require a method that “should use only information in the matrix of 
the system and as little extra information as possible [Ref. 11].” In other words, 
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we want to extend all of the benefits of multigrid methods to a new type of solver, 
algebraic multigrid. 
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III. ALGEBRAIC MULTIGRID 



Algebraic multigrid (AMG) methods are a required extension of the multigrid 
methodology. They were first introduced in the mid 1980’s by A. Brandt [Ref. 12], S. 
McCormick and J. Ruge [Ref. 13] and further developed by J. Ruge and K. Stiiben 
in their paper, Algebraic Multigrid [Ref. 14]. The methods arose from the need to 
solve purely algebraic systems ( matrix equations) of the form 

n 

Au = f or Y a ij u j = fi (i = 1, . . . , n) (III.l) 

j = i 

with the “multigrid-style” cycling process without any underlying geometry. Al- 
though there is nothing special about the form in (III.l) as is, it is the form required 
for our work in AMG. In their paper, Ruge and Stiiben show the usefulness of AMG 
when applied to a “model class” of matrices - symmetric M-matrices. (Recall that 
the matrix A is an M-matrix if an > 0 and a,y < 0 for i ^ j.) They also relaxed 
the assumption to weak diagonal dominance (|o tt | > K'jl, i = l,...,n) and 

showed uniform two- level convergence for that case as well. In addition, they imply 
that AMG may be attractive as a “black box” solver for algebraically posed elliptic 
problems as well as for certain other types of operators that generate matrices with 
similar characteristics. 

Moreover, for particular types of matrix equations, algebraic multigrid is a 
robust and efficient matrix equation solver. It is designed to solve these types of 
matrix equations (III.l) using the same principles of standard multigrid methods, 
without requiring the need of an underlying geometry of the continuous problem. 
Application of AMG to many problems involving systems that generate symmetric 
M-matrices is already shown to be efficient. We discuss below some other types of 
problems to which AMG can successfully be applied. 



27 



A. WHY ALGEBRAIC MULTIGRID? 



AMG, like multigrid, uses the same multi-scale hierarchy to solve the problem 
Au = /. At each level during the solution process, the system is “solved” and the 
error corrected until the process terminates at the lowest level (when no more levels 
are required) and then the system is solved by using one of a host of direct methods 
such as Gaussian elimination. Thus far, there seems to be little difference, but shortly, 
we will see that there is. 

It is a solution method that is ideal for solving more general large, sparse 
algebraic equations because, unlike MG, it is not dependent on particular domains 
or operators. Although AMG can solve standard elliptic problems on uniform rect- 
angular grids, for such problems, it is a less efficient solver than highly specialized 
geometric multigrid solvers. However, once the setup phase (discussed in the next 
section) is completed, AMG is quite competitive. But when geometric grids become 
impossible to design and MG cannot be tailored to fit the problem, this is when AMG 
should be used. It has been shown in various papers such as [Ref. 13, 14, 15] to have 
favorable results for these more complicated domains. 

Even though originally designed to solve large, sparse matrix equations that 
involved M-matrices, application of AMG to other types of not-so-sparse matrices has 
increased in recent times. This is primarily due to current problems of interest which 
are on unstructured grids and are extremely large, usually on the order of millions of 
equations and unknowns. In fact, in some particular cases, algebraic multigrid has 
been proven to be the fastest method available. One problem, for example, is the 
transistor placement problem in integrated circuit design. This particular problem 
concentrates on layout optimization and solves a system that determines the optimal 
location for hundreds of thousands of cells, which consist of hundreds of transistors, 
on the chip surface. Regler and Riide show that AMG is a competitive alternative to 
current methods in layout optimization. 
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B. APPLICATIONS OF AMG 

Algebraic multigrid methods are useful when it comes to solving certain sys- 
tems that either contain no geometric structure or that produce irregular grids tailored 
to complex domains or operators. When the application of standard multigrid meth- 
ods proves to be impossible or entails a high degree of difficulty on certain problems, 
algebraic multigrid may be a better choice. These types of problems are described in 
more detail in [Ref. 14] but here, they are summarized. In general, AMG should be 
considered over geometric multigrid 

(1) when the problem involves complex domains so that any discretization we 
choose is still too fine to serve as the coarsest grid. In this case, the work required 
in geometric multigrid to determine the interpolation schemes would be prohibitive. 
AMG is clearly the choice here since the coarsening process is automatic. In fact, for 
this case, AMG performs quite favorably in terms of operation count and CPU time 
[Ref. 15]. 

(2) when uniform coarsening is not allowed at all on the finest grid. This 
happens, for example, when a finite element discretization is applied using irregular 
triangulations. Since AMG requires no geometric structure (i.e., no uniform grids), 
it is ideal. 

(3) when the interpolation operator is chosen so that it becomes very difficult to 
find a relaxation process that complements the operator in smoothing the error enough 
to allow a decent coarse grid correction. Sometimes, it is impossible to determine a 
relaxation process at all. AMG is not affected here since interpolation is defined by 
the entries in the matrix, and the weights are chosen tot reflect the relative strength 
of each connection. 

(4) when a problem is entirely discrete, especially if the problem contains no 
geometric structure at all. Since AMG does not depend on an underlying continuous 
problem, it is clearly the method of choice. This applies to problems where we are 
only given the system matrix A and the right-hand side vector /. 
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C. THE STEPS OF AMG 

Algebraic multigrid methods are designed with regular (geometric) multigrid 
methods in mind, in that they use the principles of multigrid but do not require nor 
rely on the geometry of the particular problem being solved. Instead, they explicitly 
use information from the matrix of the system, i.e., the coefficients of the unknowns. 
Application of AMG to the problem Au = / requires a two-part process. There is 
an initial step which defines the intergrid transfer and coarse-grid operators and, in 
addition, chooses the coarse grid itself. This first step of the process is called the 
setup phase and most of the work involved in writing the algorithm takes place in 
this portion of AMG. “The generality of AMG must be paid for by a setup phase 
that may take 80% or more of the overall time [Ref. 15]”. The second part of the 
AMG process is called the solution phase. In this part, regular multigrid cycling is 
performed on the components of the matrix until a predetermined tolerance is met. 

1. The Setup Phase 

As defined in [Ref. 14], a brief outline of the algorithm used in the setup phase 
on the problem Au = f is provided below: 

1. Set m = 1. 

2. Choose the coarse grid and define /™ +1 , the interpolation operator. 

3. Set /;+> = (/” +1 ) T and A"+> = C + M”/” +1 . 

4. When fl m+1 is small enough, set q = m + 1 and stop. Otherwise, let 

m = m + 1 and return to step 2. 

Here, we define q to be the number that represents the coarsest grid. This number is 
used in the theoretical bound on the asymptotic (V-cycle) convergence factor, p. In 
practice though, p is ^-independent and by increasing the accuracy of interpolation 
for smooth errors, we find that AMG is quite an efficient method. Therefore, the large 
investment in writing the setup phase may save time (and cost) in the setup work for 
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later problems that use the same, or nearly the same, matrix. This effectively reduces 
the time required for the solution of follow-on problems. 

2. The Solution Phase 

The solution phase is part of the algorithm in which regular multigrid cycling 
takes place. It involves a relaxation process to remove the oscillatory modes of the 
error and an intergrid transfer operators to move the problem to different levels. This 
phase of AMG, when the operators and relaxation are chosen properly, is a very 
efficient solver for the problem on the finest grid [Ref. 14]. 

D. DEFINING THE “GRID” 

One of the main ideas of algebraic multigrid is to split the set of variables 
on a particular level into two disjoint subsets. The first set contains coarse level 
variables (C-variables) and the second is the complementary set of fine level variables 
(F- variables). The fact that we require the sets to be disjoint will be important as we 
will see shortly. If we now define h and H to be any two consecutive levels with h as 
the fine one then we can set C h — {j \ j € C} and F h — {j | j Gf}. Furthermore, 
let i £ fl h = C h U F h be defined as a point in a fictitious plane. In this sense, i 
is nothing more than a reference to the variable uf. Now, for example, A h u h = b h 
can be interpreted as a grid equation on a fictitious grid Q, h . Furthermore, we define 
Q, h _ ph an( ^ qH _ qH ^he d e t a ii i n defining the grid can be quite cumbersome 
and it is not our intention here to reproduce what can be found in [Ref. 14]. 

E. CONNECTIONS AND CONVERGENCE 

Next we introduce a new term called connections. Connections refer to the 
relationship between grid points in the sense of a directed graph which is associated 
with any matrix. On any level h , a point i £ £l h is defined to be (directly) connected 
to a point j £ £l h if a,ij ^ 0. Furthermore, we define the (direct) neighborhood of a 
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point i G tl h by 

N ?={ j € n K \ (111*2) 

Let us now concern ourselves with interpolation along direct connections. In 
this case, we consider only the interpolation points C; with Ci C N{ D C and the 
corresponding weights 



Wik = Vi 




(i € F, k£ Ci) 



(111*3) 



with 



0 < rji < 1 / E K • 

leCi 



(HI-4) 



In order to maintain two-level convergence, it is sufficient to require for every i € 
F, k eCi 



where s,- < 1, which is ensured by (III.4). From (III.5), we can easily derive more 
practical conditions to effectively develop an automatic coarsening algorithm that has 
/5 as an input parameter. The algorithm will also choose interpolation with weights 
(III.3) satisfying (III. 5). 

In order to understand the effect of /3 in (III. 5), we state the following result 
presented in [Ref. 14]. 



Theorem: Let [3 > 1 be fixed. Assume for any symmetric, weakly diagonally 
dominant M-matrix A h the C-points are picked such that, for each i E F, 
there is a non-empty set C,- C A,- fl C with 

fi £ < > 4 (me) 

itc, 



and define the interpolation weights (III. 3) by r/,- = 1/ a ij ■ Then two-level 

convergence is satisfied. 



For a given /?, condition (III. 6) is satisfied by many choices of the sets C 
and Ci ; but regardless of how we chose those sets, we always want them as small as 
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possible. This condition requires that we make each F-point i strongly connected to 
its interpolation points. In other words, we actually want each |a,j| (j (E Ci) to be 
comparable in size to the largest of the |a t -*| ( k € iV,) or equivalently, 

\a;A « 7 max la,*! for j 6 C ; and 0 < 7 < 1 . 

J keNi —I — 

The assumption in (III. 6 ) gets weaker as /3 gets larger; but, for large /3, (III. 6 ) 
allows for rapid coarsening. The price of this though is slower two-level convergence 
speed. Although, when /? = 1, we achieve the best convergence but the expense of 
the method becomes extreme. Ruge and Stiiben explain that the best compromise 
is when (3 = 2. This is when about 50% of the total strength of connections of every 
F-point will be represented on the coarser level. 
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IV. CONSTRUCTING THE 
INTERPOLATION OPERATOR 



The choosing of an interpolation operator is “dictated by the desire to achieve 
good theoretical estimates of the convergence of the iterations as well as by bounds 
on the computational complexity of the algorithm [Ref. 11].” It is of the essence 
in defining the interpolation operator that the range of interpolation approximates 
those errors not effectively reduced by relaxation. To make it successful then, AMG 
must be constructed in such a way to ensure that this automatically happens. 

A. ALGEBRAIC SMOOTHNESS 

For the problem Au = /, let us define a smoothing process 

ni m . rm ( a m\ — 1 rm 

u new = Gr U old + / — G {A ) f 

where G is a (linear) smoothing operator such that G m : 1Z nm — ► 1Z nm . In AMG, an 
error e is said to be smooth if HGeHj « He]^. In other words, the error is smooth when 
it is slow to converge with respect to the smoothing operator. If we now assume that 
e is a smooth error then the residual r = Ae is small compared to e. This occurs in 
most of the common relaxation processes. Knowing this, we expect |r,| <C an |e t -| for 
all i and so for (algebraically) smooth error we have Ae ~ 0. A good approximation 
for each e, can then be obtained from 

Ae = anei + ^ a ij e j ~ 0 

jeNi 

or 

ei = — a v e J ( IV -1) 

a ” it* 

where the neighborhood AT,- of point i is defined in (III. 2). 

With this fact, it is easy to construct an interpolation operator which guaran- 
tees that the smooth error lies in the range of interpolation. One standard operator 
might be 
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0 H 



E Wik^k 



H 



l k€C h 

with Wik = 6ik if i £ C h . Recall that 



for i G C = Sl H 
for i G F = n h 



£>ik = 



0 

< 

1 



for i =£ k 
for i = k. 



However, this basic “model” interpolation operator isn’t very effective since it 
is not usually “local” enough. Here we define “local” to be wik = 0 for all i G F h 
unless k € Cf, where C- 1 C C h . For reasons of efficiency, we require that Cf be 
some small set of interpolation points. When working with symmetric M-matrices, 
[Ref. 14] shows for the case when l a »jl ~ a a that smooth error generally varies 
slowly in the direction of strong connections. This fact will be important later on 
when constructing weight definitions i 0^ for the interpolation operator along those 
connections. 



B. INTERPOLATION ALONG DIRECT CONNECTIONS 

Let us consider only those interpolation points C,- such that C, C Ni fl C. 
The interpolation weights then have the form a\ k | for i € F and k € C; 

where 0 < 77 , ■ < (j2ieCi |°^|) • require this definition to ensure that wu- < 1 . 

Let us take a look at the case when the matrix A is strictly an M-matrix (with no 
assumption on symmetry), i.e., an > 0 and a,j < 0 for i ^ j. If we set Ci = Ni, then 
(IV. 1 ) gives 






*13 | 



i€iV, a <« 



e j' 



(IV.2) 



Now we choose wn- = a-^j /an. However, this would imply that C = Cl but we want 
Ci to be small; so generally, we desire the set Di = Ni — C,- ^ 0. What we want to do 
is to “distribute” the non-interpolatory connections a t y ( j G Di) to the interpolatory 
point an- 
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If \a,ij\ is small, then adding a,j, for j e A, to the diagonal element an should 
not hurt accuracy very much since the contribution from a,j is negligible. If, however, 
|a,j| is large, then we know this replacement is reasonable since smooth error generally 
varies slowly in the direction of strong connections. To determine the error in this 
case, we begin with (IV. 1). 



then 



or 



an&i — ^ 1 aij Cj 

j€N, 

— — aijtj — 'y' GijGj 

jec, jeDi 

s ■■ 1 V ' 

set ej»e t 

~ dijtj — a ij e i 

jeCi jeD t 



I an + ^ ^ j ~ ^ y CLijCj 

\ j€D t / jeCi 



Zi = 



jeCi 

jeDi 

E lay \ej 
jec, 

“I" y . aij 

jeDi 



(IV.3) 



Since A is assumed to be an M-matrix, the denominator in (IV.3) is positive. We 
now have a good representation of the error. 



C. THE ITERATIVE WEIGHT DEFINITION 

Let us now consider the entries in an arbitrary matrix A. Figure 12 shows a 
graphical relationship between various points in the matrix. We will use this illus- 
tration in our derivation of the weights. Recall our earlier assumption that Ae « 0 
when e is smooth. Then, for a smooth error and for a point i € F, we have 

anti = - Y Q ikek ~ Y a d e i “ Y aiie ‘ ( IV - 4 ) 

keCi jeF ieN„ lec, i$c. 
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le Nj 




Figure 12. Relationship of fine grid to coarse grid points in the matrix. Points i and 
j are the fine grid points; the others are coarse grid points, both weak and strongly 
connected, interpolatory and non-interpolatory. 

and similarly for a point j € F , we have 

®jj^j = ^ ^ ^ ^ ajiCi ^ ^ (IV. 5) 

ArGC'j i£F m£Nj , m£C, m$.C 3 

To make the interpolation operator “small,” we must distribute the non-interpolatory 
connections (a,j for j ^ C,) to the interpolatory points. In basic terms, we need to 
approximate ej ( j 6 F) and e; (/ 6 N{, l € C, l £ Cj ) in terms of ( k € Cf). 
Likewise, we need to approximate e, ( i € F) and e m (m € Nj, m £ C, m ^ Cj) in 
terms of ek (k 6 Cj). 

Let us assume that the interpolation weights have been determined for the 
interpolatory points k € Cj. That is, the interpolatory weights to define ej for 
j € F are known. Then ej could be approximated by a weighted average (for use in 
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approximating e;): 



e j 



Wjkl^kl + Wjk 2^2 
w jki "f" Wjk 2 



so in general 



trj ^ 



E w jk e k 
kec, 

E Wjk 

k€C, 



(IV.6) 



Then for i £ F, we substitute tj in (IV.6) into ej in (IV. 4). This leads to 



diiCi — djk^k djj 

k€Ci jeF 



^ E Wjk&k 



kec. 



E wj k 

\ keCi 



X o-ii^h 
ieN it tec, itCi 



(IV.7) 



while for j € F, we replace e t - in (IV. 5) by the equivalent approximation for e, in 
(IV.6). This gives us 

( E Wik& k \ 

a jj£j = X! a jk&k X! 



AreCj 



t'eF 



kec j 

E 



E 

m£Nj, m£C, mfcC 3 






(IV.8) 



If we look at the two previous equations closely, we see that (IV.7) is being 
used to define W{ k which are used in (IV.8) and (IV.8) is being used to define Wj k 
which are used in (IV.7). Thus, the interpolation weights at point i are dependent 
on the interpolation weights at point j, and vice-versa. Hence, we have an implicit 
system to define the interpolatory weights. This implicit system governs how the 
weights are computed and is the basis for the iterative weight definition (IWD) of 
AMG. 



D. INTERPOLATION CONSTRUCTION USING THE 
ITERATIVE WEIGHT DEFINITION 

This section presents an algorithm for coding the iterative weight definition of 
algebraic multigrid. We begin by starting at a fine grid point, say point i, and dividing 
its neighborhood, N{, into four distinct groups. Figure 13 represents an illustration 
of a typical neighborhood of an arbitrary point i. It also shows the elements in each 
of the groups to which we will refer in this discussion of construction. 
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We begin the algorithm by dividing the entries of the matrix into four groups. 
For each point i, divide N{ into four groups: 

le N; 




Figure 13. Illustration of the neighborhood (N{) of a fine grid point i and its relation- 
ship to other fine and coarse-grid points in the matrix. For each point i, N{ is divided 
into four distinct groups. 



Group A 



Group B 



Group C 



Group D 



Coarse-grid interpolatory points, C,. In Figure 13, these corres- 
pond to points &i, k 2 and k 3 . 

Coarse-grid and fine-grid non-interpolatory points that are weak- 
ly connected to point i. In Figure 13, these are points m and l. 
Fine-grid points strongly connected to point i. In Figure 13, this 
is point j. 

Coarse-grid points not weakly connected to point i but not in 
Ci . In Figure 13, this is point n. 
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1. Initialization 

Step 1: Initialize the variables (vectors): 
For all i E F do 
set 

Ci = Si DC 



set 



define 



end 



w ik = 



&ik 



Y a il 

ieCi 



Cli = {j | dij is small } 



Here, we note that Si is the set of points strongly connected to point i and fl, is the 
set of points weakly connected to i. 

2. Calculation 

Step 2: Calculate the interpolatory weights, to,*, by using the approximation 
At « 0: 

For each i € F do 

anti ft: 





(IV.9) 


keCi 




aiiei 


(IV.10) 


l€Qi 




yi a tj e j 


(IV.11) 


jeF, j$n t 




^ ] ^in^n 


(IV. 12) 


n£C, n£Ci, n£Q, 





approximate the errors: 

for group B points, distribute the quantity in (IV. 10) to the 
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diagonal by the approximation 



ei « e,-. 

for group C points, approximate 



(IV. 13) 



e . /■'o 

3 



E w ik e k 

kec, 

E w ik 

keCi 



(IV. 14) 



but distribute the quantity in (IV. 11) to the diagonal only 
when Wi k > 0. 

for group D points, approximate 



E Gik^k 

kec, 

E ® ik 

kec. 



(IV. 15) 



but distribute quantity in (IV. 12) to the diagonal only when 
a,j is of the right sign. 

Either stop here, or 
For each i € F do 
set 



Ci = {j | Wij > 0} 

goto Step 2. 

end 

end 

We give some insight to the iterative nature of this last step. The routine ends 
when there are no more fine-grid points or the set of interpolatory coarse-grid points, 
Cj, is empty. If this is not the case, we check the next fine-grid point i and refine the 
set Ci to include only those weights that are of the correct (positive) sign. If all the 
weights are now of the correct sign, then we end, if not we refine Ci again and repeat 
the loop until we meet an end of loop criteria. 
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E. INTERPRETATION OF THE ALGORITHM 



In this section, we interpret the code of Section D. But before we begin, we 
make a few simplifications. Define the sets 

A = { k\k€Ci } 

B = {/|/ea> 

C = {j\j € F,j i ft,-} 

D = {n|rc e C, n £ Ci,n £ OJ . 

Now, the equation involving (IV.9)-(IV.12) becomes, 

dii&i = dj k G k d" ^ ' djj&l -J- 'y * QjjGj -)- ^ ^ dj n e n j . (IV. 16) 

\ A B C D / 

In the code, the sum over the set B becomes after substitution of (IV. 13) into (IV. 10) 



(IV. 17) 

B 

the sum over the set C becomes after substitution of (IV. 14) into (IV.ll) assuming 

w ik > 0 



E w ik e k \ 

A (IV. 18) 

and the sum over the set D becomes after substitution of (IV. 15) into (IV. 12) assum- 
ing a,j is of the right sign 




After interchanging sums in 



'y ^ a in 

D 



E a i k e k \ 

A 





(IV. 18), we get 



£ 



Edij \ 

C 



E Wik 

k£C t 



Wik&k 



(IV. 19) 



(IV. 20) 
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and after switching the sums in (IV. 19), we get 



E 



E«m \ 

_D 

E a ik , 

\ kec, } 



O'ik^k* 



(IV. 21) 



With (IV. 17), (IV. 20) and (IV.21), (IV.16) becomes 



^ ^ CLjiCj ~ 
B 

or equivalently 



^ ^ ajk^k ^ ^ 
A A 



( £ a ij \ 

C 



W. 



‘k^k H" ^ ^ 



/ Eflin \ 



D 



£ O'ik 
\keCi ) 



CLik^k 



| Q>ii ^ ^ Uj/ I C{ ~ ^ ^ 

V B ) A 

We now rewrite (IV. 22) as 



aik + 



E a ij 

C 

, E toi* 

\keCi 



Wik + 



E G m \ 

J} 

£ aik 

\kec t J 



a{k 



e k . (IV. 22) 



e i ~ £ 

,4 



£ a *i 

a * fc + l 

k<=C t 



E a - 

Wt'fc + ( I a ik 

Kk GC t - 



(flit + E 

V B 



<2j; 



e fc . 



The term in the braces is W{ k . So our formal definition of the iterative weight definition 
is 



aik 



£ a *j 

c 



TT<' 



nOld 



old w ik 



W7L - 



in new ~ 
W ik — 






£a, n 

D 

t 

^fc€C t 



aik 



( a ti + E a i/) 
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V 



NUMERICAL RESULTS 



In this chapter we present some numerical results on problems involving dis- 
cretized partial differential equations as well as non-geometrically generated matrix 
equations using algebraic multigrid methods. As will be seen shortly, AMG methods 
using the iterative weight definition (IWD) may not be robust but are more efficient 
on certain problems than standard AMG weight definitions. 

Before we begin with the results, a few definitions are in order. We define a 
strong connection as 



|a«il = 7 max |a ifc | for j € C; 

k£N t 

where 7 is a parameter that varies: 0 < 7 < 1. The choice we use in our experiments 
to define a strong connection is 7 = 0.25 which has been experimentally shown to 
yield good results. Most problems presented use this choice of 7; however, for some 
problems we test with 7 = 0.5 because, in certain cases, favorable results have been 
found (and is another avenue of research). For those problems, it will be specifically 
stated what 7 is. For reference purposes, (1,1) V-cycles with Gauss-Seidel relaxation 
and C/F-ordering of points were used in all tests. The convergence factor ( p ) was 
computed by 



fail 

vim; 



(V.l) 



where n equals the number of cycles performed, Rf is the final residual and is the 
initial residual. The following notation will also be used throughout this chapter: 



a A Operator complexity 
a n Grid complexity 

p the asymptotic convergence factor of the 20 th V-cycle 



where a A is the ratio of the total number of nonzero entries in all the matrices to that 
of the number of nonzero entries in the fine-grid matrix and is the ratio of the 
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total number of points on all grids to that of the number of points on the fine-grid. 
o A and are other measures that describe the performance of AMG. After the first 
few experiments we use p instead of p since the 20 tA V-cycle is more representative 
of the true convergence rate. 



A. ISOTROPIC PROBLEMS 

We wish to establish some baseline behaviors of each method in the absence 
of any irregularities. We discretize the Laplace operator (V 2 ) on the unit square with 
Dirichlet boundary conditions using several different finite difference discretizations. 
For each problem, we use a uniform grid with mesh size h = 1/64. If we let N = 1/h 
then we create a square matrix of size N 2 x N 2 (or 4096 x 4096 for our problems). 
The stencils we use are as follows: 



( 1 ) 



h 2 



-1 

-1 4 

-1 





> 



( 3 ) 



1 

Sh 2 



-1 


-1 


-1 


-1 


oo 


-1 


-1 


-1 


-1 



( 4 ) 



20h 2 



-1 


-4 


-1 


-4 


20 


-4 


-1 


-4 


-1 



These stencils arise naturally from the discretization of 

— V 2 u(z,j i) = f(x,y) in ft 
u = 0 on dS l. 

To get stencil (1), for example, we replace the derivatives of 

U xx Uyy — f(x,y) 



(V.2) 



(V.3) 



by second order finite differences. If we let h x = 1/N and h y = 1/M , where N and 
M are positive integers then 



^ XX 



— («,•_ i,j - 2 mj + u i+1J ) 



(V.4) 
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and 



u 



yy 



7 9 ^UiJ + Uij+l) . 






(V.5) 



Substituting (V.4) and (V.5) back into (V.3) and with the definition /,j = /(;r,-,t/j), 
then (V.2) becomes 

+ ^i+lj 2 ti;j + Uj‘j+1 

— Ji,j 

(V.6) 






*2 



U 0 j ^N,j — M«,0 — — 0 



where 1 < z < iV — 1 and 1 < j < M — 1. 
In terms of a stencil then, we get 



- V 2 = 



-1 

-1 _ 2 _ _ 2 _ -1 

A? + If 

-1 



If we now let h = h x = h y , then (V.7) becomes 



-V 2 = — 

h 2 



-1 
-1 4 

-1 



-1 



(V.7) 



which is exactly what we were trying to show. The other stencils can be similarly 
derived. 

Table I compares the results of the iterative weight definition against standard 
AMG weights on Laplace’s equation (1.1) using the stencils above. In all tests the ini- 
tial guess for u(x,y) was randomly generated. On Laplace’s equation using standard 
stencils discretized using finite differences, we note that there is little, if any, differ- 
ence in p of the iterative method over the standard definition. As expected though, 
AMG produces excellent results on these model problems using either method. 
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Table I. Convergence rates in solving Laplace’s equation using different isotropic op- 
erators discretized using the finite difference method. 





Standard weights 


Iterative weights 


Stencil 


P 


Complexity 


P 


Complexity 


<7^ 


<7 U 


a A 


<r“ 


1 


0.05037 


2.3714 


1.7009 


0.05037 


2.3726 


1.6990 


2 


0.07903 


2.2475 


1.6799 


0.07903 


2.2588 


1.6797 


3 


0.08409 


1.4162 


1.3464 


0.08109 


1.4172 


1.3440 


4 


0.06541 


3.7120 


1.9963 


0.07938 


3.9734 


2.0522 



B. THE ANISOTROPIC PROBLEM USING THE FI- 
NITE DIFFERENCE METHOD 

Now that we have some baseline results of both methods in the absence of any 
irregularities, we wish to try a new problem, 

—eu xx — u yy = 0 in Cl 
u = 0 on dLl 

where the operator, —eu xx — u yy , is the anisotropic Laplacian. Results in this section 
are given for various levels of anisotropy (e) and again we use second order finite 
differences to discretize the problem. This section will demonstrate the ability of 
AMG to tailor the coarsening algorithm to the individual problem. Other results can 
be found in [Ref. 9]. The domain in the following problem is again the unit square 
with Dirichlet boundary conditions. The problem is discretized on a uniform grid 
with h = 1/64, using the 5-point stencil 




h? 



-1 

— s 2 T 2s — s 
-1 



(V.9) 



We will now apply both weight methods to problem (V.8) with f(x,y) = 0. 
The motivation behind this test is to have knowledge of the error. With the error 
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known exactly, it is easier to compare the true performance of both methods. As 
before, the problem is discretized on a uniform grid with h = 1/64, using the 5-point 
stencil in (V.9) above. This generates a 4096 x 4096 matrix. The convergence rate in 
Table II is that of the 20th iteration (/>) vice the average defined in (V.l). In Table 
II, we see that the iterative weight definition produces a better convergence rate over 
standard AMG weight methods for all values of e in the anisotropic problem but these 
improvements are utterly insignificant since they are so small. In addition, IWD does 
not always produce better results in terms of solver complexity or grid complexity. 
Clearly the standard weight definition is the better in this case since the increased 
cost associated with IWD makes it unfavorable when the results are practically the 
same. 

Table II. The anisotropic problem of the Laplacian operator discretized using the finite 
difference method. 





Standard weights 


Iterative weights 


6 


P 


Complexity 


P 


Complexity 


a A 




a* 




0.0010 


0.04549 


2.6319 


1.9644 


0.03944 


2.6299 


1.9609 


0.0100 


0.06731 


2.7707 


1.9590 


0.06698 


2.7707 


1.9590 


0.1000 


0.05674 


3.3067 


1.8655 


0.05672 


1.3103 


1.8655 


0.5000 


0.06349 


2.4027 


1.7068 


0.06085 


2.3272 


1.6926 


1.0000 


0.05721 


2.3714 


1.7010 


0.05037 


2.3727 


1.6990 


2.0000 


0.06951 


2.4271 


1.7129 


0.05206 


2.3338 


1.6956 


10.000 


0.05576 


3.5391 


1.9033 


0.05575 


3.5691 


1.9070 


100.00 


0.06753 


2.8047 


1.9658 


0.06720 


2.8043 


1.9653 


1000.0 


0.04544 


2.6319 


1.9644 


0.03946 


2.6299 


1.9609 



Again, one can see in Table II that the differences, if any, between the methods are 
minimal. For these simple problems presented thus far, we see little benefit in using 
the iterative weight definition. In fact, if we consider the cost associated with IWD, 
we conclude that the method is worse since the standard weight method is a less 
expensive routine. 
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C. THE ANISOTROPIC PROBLEM USING THE FI- 
NITE ELEMENT METHOD 

We present results in this section using the discretized Laplacian operator. 
Whereas Sections A and B involved discretizations using the finite difference method, 
we apply the finite element method to see if there are any improvements. 

Table III provides results generated by elongating the discretized grid for 
Laplacian operator in the x-direction. We also compare the convergence rates of 
both methods as we change the definition of strong connection from 7 = 0.25 to 
7 = 0.5. In the table, dx = edy means that for every point in the y-direction (dy), 
there corresponds e more points in the x-direction {dx). The intention of this test 
is to see how well the iterative method tailors to semi-coarsening to improve the 
convergence rate. 



Table III. The anisotropic problem with zero RHS discretized using the finite element 
method. 





standard weight 


iterative weight 


size {N x N ) 


dx — edy 


7 = 0.25 


IO 

O 

II 


7 = 0.25 


IO 

O 

II 


matrix 


grid 


6 


P 


400 


20 


25 


0.14 


0.14 


0.14 


0.14 


900 


30 


10 


0.45 


0.23 


0.13 


0.13 


900 


30 


100 


0.83 


0.53 


0.82 


0.23 


1225 


35 


50 


0.25 


0.14 


0.15 


0.14 


10000 


100 


10 


0.47 


0.24 


0.14 


0.14 


10000 


100 


100 


0.93 


0.55 


0.93 


0.28 



From inspection of Table III, it is clear that for certain problems, the iterative 
weight definition is significantly better. In some cases, IWD is better by over | and 
the extra cost associated with using this method has more than paid for itself with the 
large reduction in convergence rate. The natural questions then are why is it better 
on some and not on the others? What are the special properties of the problems in 
which IWD outperforms the standard weight definition? Does semi- coarsening play 
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a role? Before we can answer these questions, we look at some other problems with 
domains that are much more complicated. 

D. PROBLEMS WITH COMPLEX DOMAINS 

In this section, we present some results for some complicated domain problems 
using the standard Laplacian operator 

~ dx 2 + dx 2 ' 

As in Sections C, we use the finite element method to discretize the problem. Figures 
14, 15 16, 17, 18 and 19 are graphic illustrations (meshes) of the problems we inves- 
tigate. The size of each problem varies and is noted in the following tables. Table 
IV shows the results of some initial experiments on these more complicated domains. 
Much like the problems we have already seen, we expect that little will be gained by 
using the iterative weight definition since on standard Laplacian operators, regular 
AMG already works great. 

1. The Meshes 




Figure 14. Mesh 1. 
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Figure 15. Mesh 2. 




Figure 16. Mesh 3. 
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Figure 17. Mesh 4- 




Figure 18. Mesh 5. 
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Figure 19. Mesh 6. 



2. Results of Complex Domain Problems 



Table IV. Convergence rates of problems high in complexity using the standard Lapla- 
cian operator discretized using the finite element method. 





Standard weights 


Iterative weights 


Mesh 


Size 


P 


Complexity 


P 


Complexity 


a A 




<7^ 




1 


561 


0.8717 


2.1990 


1.7184 


0.8496 


2.1998 


1.7184 


2 


513 


0.0896 


2.1250 


1.6530 


0.0902 


2.1118 


1.6472 


3 


321 


0.0766 


3.1268 


1.9595 


0.0746 


3.0207 


1.9315 


4 


757 


0.1386 


2.9944 


1.9511 


0.1779 


2.9780 


1.9406 


5 


1256 


0.1460 


3.3765 


1.9936 


0.1460 


3.3000 


1.9745 


6 


1080 


0.1501 


2.3064 


1.7630 


0.1399 


2.3541 


1.7769 



With these results, we suppose that with simple operators such as the standard 
Laplacian or even the slightly skewed Laplacian the iterative weight definition will 
not perform better than regular AMG methods. This is largely due to the fact that 
conventional MG works extremely well on standard (elliptic) operators. We must test 
our new method using a more irregular type of operator to see if IWD will fare better 
than standard weight definitions. 
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E. COMPLEX DOMAINS AND OPERATORS 



Laplacian operators have not seemed to produce results in favor of the iterative 
method so we hypothesize that IWD may outperform the standard weight definition 
method when the operator is more complex. To get a better idea of what types of 
operators or on what types of problems the iterative weight definition improves p. 
we develop a more complex operator and reevaluate some of the same problems as in 
Section D. 

Specifically, we experiment with those geometries as illustrated in Figures 14, 
15, 17, and 16. The size of each problem here also varies and may be different from 
those presented in Table IV. The sizes of each problem are included in the tables. 
For problems in this section, we use the irregular operator 



L = * x+y ir- 2 + (* + V): 



+ e x+y -^- + i-. 



dx 2 dy 2 ' dx ' dy 

Table V provides results of both methods using the definition of strong con- 
nection 7 = 0.25 whereas Table VI uses 7 = 0.5 as the definition. 



Table V. Convergence rates of complex problems using an irregular operator discretized 
using the finite element method with 7 = 0.25. 





Standard weights 


Iterative weights 


Mesh 


Size 


P 


Complexity 


P 


Complexity 




<7“ 


<7" 


u* 


1 


561 


0.1367 


2.3281 


1.7558 


0.1432 


2.3891 


1.7825 


2 


513 


0.0984 


2.3862 


1.9415 


0.1279 


2.3543 


1.9376 


3 


1249 


0.5912 


3.6256 


2.0200 


0.2877 


3.6133 


2.0168 


4 


2889 


0.7949 


3.4115 


2.0048 


0.8880 


3.2519 


1.9740 



Tables V and VI and reveals something interesting. The iterative weight defi- 
nition significantly improves the convergence rate over the standard weight definition 
by approximately |, but only in one case. What is special about that case? Clearly, 
it is not this irregular operator alone nor the specific problem but a combination of 
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Table VI. Convergence rates of complex problems using an irregular operator dis- 
cretized using the finite element method with 7 = 0.5. 





Standard weights 


Iterative weights 


Problem 


Size 


P 


Complexity 


P 


Complexity 


< 7 ^ 


a u 


a A 


< 7 “ 


1 


561 


0.2700 


2.7205 


1.9287 


0.1194 


2.7501 


1.9305 


2 


513 


0.1106 


2.4271 


1.9552 


0.2863 


2.4334 


1.9591 


3 


1249 


0.7573 


3.0000 


1.9976 


0.4668 


3.0286 


2.0120 


4 


2889 


0.9013 


2.9619 


2.0142 


0.9400 


2.9612 


2.0059 



both that holds the answers to our quest for classification. Moreover, in the other 
three problems, IWD performs worse. 

Let us now investigate the case (mesh 3) where IWD works the best to see if 
there are any clues to its usage. Figure 20 shows not only the sparsity pattern of the 
matrix generated but more to the point, the locations of the entries in the matrix. 




Figure 20. The sparsity pattern of mesh 3. 

Figure 21 shows where the positive and negative entries are. Together, these 
figures reveal that mesh 3 is not a M-matrix and so we have a non-M-matrix case where 
the iterative weight definition should be used over regular AMG weight definitions. 
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Figure 21. Location of positive and negative entries of mesh 3. 



Table VII shows the difference between the iterative weight definition and the 
standard weight method over a sequence of refinements on that particular problem 
(mesh 3). In all cases, we see a major improvement in convergence rate. 



Table VII. Convergence rates of the problem in Figure 16 using an irregular operator 
discretized using the finite element method. Different refinements to the mesh are 
presented. 





Standard weights 


Iterative weights 


Refine 


Size 


P 


Complexity 


P 


Complexity 


<r A 




a A 


<x“ 


1 


85 


0.3967 


2.6373 


1.9176 


0.3411 


2.7199 


1.9412 


2 


321 


0.4368 


3.3284 


2.0125 


0.2684 


3.4382 


2.0467 


3 


1249 


0.5912 


3.6256 


2.0200 


0.2877 


3.6133 


2.0168 


4 


4929 


0.7351 


3.7220 


2.0016 


0.3582 


3.6430 


1.9846 



Table VIII shows other experiment involving this irregular operator on the 
problem (mesh 5) shown in Figure 18. 
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Table VIII. Convergence rates of the problem in Figure 18 using an irregular operator 
discretized using the finite element method. Different refinements to the mesh are 
presented. 





Standard weights 


Iterative weights 


Refine 


Size 


9 


Complexity 


9 


Complexity 




< 7 *‘ 


< 7 * 


a" 


0 


330 


1.0000 


2.9523 


1.9940 


0.9906 


2.9927 


2.0121 


1 


1256 


0.9888 


3.3843 


2.0677 


0.9881 


3.2193 


2.0239 


2 


4896 


0.8459 


3.6010 


2.0666 


0.8418 


3.4574 


2.0345 



In this case, both methods have seemed to fail all together. We conclude here that 
the operator alone isn’t the cause for the improvements seen in Table VII in using 
IWD. It must be a combination of both the geometry and the operator. 

F. CONCLUDING REMARKS 

The accuracy and efficiency of AMG in solving systems involving symmetric 
M-matrices has been shown in many papers. The use of AMG for solving other types 
of matrix equations is a topic of extensive research and variations in computational 
methodology during the interpolation phase have resulted in differing convergence 
rates. We have found that regular AMG interpolation weight definitions are inade- 
quate for solving certain discretized systems that do not lead to M-matrices. For these 
types of problems, an iterative approach to determining the interpolatory weights was 
useful. 

In applying the iterative interpolatory weight definition of AMG, we have care- 
fully balanced the requirement of keeping the interpolatory set of points “small” in 
order to reduce solver complexity while at the same time, maintaining accurate inter- 
polation to correctly represent the coarse-grid function on the fine grid. In addition, 
the extra work involved in using IWD does not significantly add to setup phase costs. 

Experimental results have shown that the iterative weight definition does not 
significantly improve the convergence rate over standard AMG when applied to M- 
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matrices, which we anticipated. However, the improvement becomes significant when 
solving certain types of complicated, non-standard algebraic equations although it is 
unclear at this stage of development what details are required to cause the iterative 
weight definition to outperform the regular weight definitions. 

We have seen that IWD is not robust; however, there are specific problems on 
which IWD should be used. When using a standard operator such as the Laplacian, 
regular AMG should be used since its performance is superb and it is less expensive 
to use than the iterative weight definition we presented here. However, when the 
operator becomes irregular and the domain more complex, IWD has been shown to 
dramatically improve the convergence rate over current AMG weight definitions and 
should be considered a viable option in the future. 
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VI. 



FUTURE RESEARCH 



Although we have shown that the iterative weight definition (IWD) of AMG 
isn’t robust, we have found certain problems where it has improved the convergence 
rate significantly. We are currently searching for more data that will enable us to 
classify the types of problems on which IWD works best. There are still many avenues 
of research and IWD is just one such direction. Algebraic multigrid research is open 
to a plethora of new and interesting ideas. The need for more efficient and robust 
solvers is never-ending and is the driving force that brought about AMG in the first 
place. To discuss all the new possibilities of improved AMG would be a paper in 
itself; therefore, we restrict our attention to that of just a few new methods on the 
forefront of improved interpolation. 

A. INTERPOLATION USING EIGENVECTORS 

One method that can improve interpolation is based on the idea of eigenvec- 
tor approximation and was originally presented in 1985 by Ruge and Stuben in [Ref. 
14]. Approximation using eigenvectors is a method used to modify interpolation and 
the coarse grid matrices on all levels. The basis of the idea is that the eigenvectors 
corresponding to the smallest eigenvalues are, in general, the algebraically smoothest 
functions and so must be better approximated by the range of interpolation. Unfor- 
tunately, the eigenvector approach requires the computation of the smallest eigenval- 
ues and their corresponding (approximate) eigenvectors. However, there are efficient 
methods for finding them and so it may not be as computationally expensive as one 
would think. Computational accuracy in determining them is rarely needed since 
linear interpolation (under the assumption that smooth functions are locally linear) 
“usually ensures accurate enough interpolation of the needed eigenvectors when stan- 
dard interpolation does not [Ref. 14].” 
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To make the method efficient, we must integrate the computation of the eigen- 
vector with the updating of the interpolation. Before the coarse-grid correction can be 
applied on each level, we must update the interpolation and the coarse-grid operator 
with the current eigenvector approximation. This must happen before the operation 
can continue to the next coarser grid. Ruge and Stiiben note that this can be done 
efficiently and that several eigenvector approximations can be calculated simultane- 
ously, if required. The concept of the eigenvector approximation is a tool used only 
for the improvement of interpolation, after which we use AMG in the standard way 
to solve the problem. 

B. THE LEAST-SQUARES IDEA 

Tom Manteuffel (University of Colorado at Boulder) and Van Henson (Naval 
Postgraduate School) recently introduced a new method that incorporates a least 
squares solution on the local level. Initial trials have produced quite favorable results. 
The idea behind this new method involves an singular value decomposition (SVD) of 
the matrix on a local level. Presented here is a brief outline of how it works. 

1 . Let i be a fine point. Define N\ to be the set of points that includes the 
fine point i and those points p connected to point i. Furthermore, define the set N? 
to include the set Ni and all the points q connected those points in N\ and continue 
this process (e.g., A 3 , A 4 , . . .) until there are no points remaining. 

2. Now select from the matrix A the rows corresponding to Ni. Remove all 
columns that contain only zeros. This removal of zero-columns yields a n x m matrix 
which we call S. For example, if we start with a nine-point stencil on regular grid, 
then we end up with a 9 x 25 matrix, S. 

3. Compute S — UEV, the singular-value decomposition of S. At least the 
last m — n (or 16, using our nine-point stencil) columns of V are null-space vectors of 
S. Let X be the m x (m — n) matrix (e.g., 25 x 16) comprising these m — n columns 
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of V. 



4. Let T be the m x m matrix whose rows correspond to the same points as 
do the columns of S. In the example of the nine-point stencil, T is of size 25 x 25. T 
is the matrix associated with solving the homogeneous Dirichlet problem on jV 2 . 

5. We now reorthogonalize the columns of X to be T-orthogonal, using a 
Gram-Schmidt process with respect to the T-inner product 

i - 1 

u i = x i ~ ( Tx h w k) for j = 1, 2, . . . , m - n, 

k = 1 



where 

Uj 

Wj = ■ , 

y{Tuj,Uj) 

and define W be the m x (m — n) matrix whose columns are u>i, te 2 , w 3t . . . , w m - n . 

6. Select I, a set of k points to be the interpolatory points from which the 
value at i is to be interpolated. Normally (but not always) the points in 1 are chosen 
from among the C-points connected to i. Permute the rows of W so that the first k 
rows contain the values of the reorthogonalized singular vectors corresponding to the 
points in I. The (fc+l) si row contains the values of the singular vectors corresponding 
to the point i. 

7. Perform k Householder reflections from the right to bring the matrix (e.g., 
25 x 16) W into the form: 



* 0 0 
* * 0 0 

* * * 0 0 

* * * * * * * 

* • • • * 



This reflection does not alter the span of 
is a unitary operation. 



<— row for interpolation point 

< — row for interpolation point 

<— row for interpolation point 

<— row coressponding to point i 

columns, since the Householder reflection 
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8. Now determine a linear combination of the first k rows that gives the first 
k entries in the ( k + l) st row. The weights of that linear combination are the desired 
interpolation weights. 

Note that the size of the entries in the ( k + l) st row beyond th e,k th column 
give a “residual” measure of how well the process works. If all the entries are “small” 
then the weights do a good job approximating the singular vectors at point i. If some 
of them are “large” then we may need to add another interpolation point to the set 
or, the set of equations selected in step 2 should be expanded to be those in jV-j vice 
those in N\ . 

Initial experiments with this method have produced excellent results on several 
standard problems. On certain problems requiring semi-coarsening (i.e., different 
coarse-grid spacings in different coordinate directions), this is the only method known 
to produce the “correct” interpolation weights. However, the method is extremely 
expensive, requiring small SVD calculations at every fine-grid point. Achieving the 
robustness of this approach at lower cost is a major focus of ongoing research. 

C. THE COMPOSITE GRID FORMULATION 

Steve McCormick (University of Colorado at Boulder) has presented some ini- 
tial work on improving interpolation by approximating “globally smooth” components 
in the neighborhood of the interpolation region. His idea is to use composite grids to 
“solve” the local problem on the coarse grids. What he recommends is to have the 
target point i and all its neighbors on the local fine grid, twice removed neighbors on 
the next coarser and more global grid, and so on, then use the current AMG coars- 
ening to solve the local problem which could in fact just be the least-squares scheme 
described in Section B. In the following discussion, L is the n x n matrix AMG is 
given to “invert” and / is the q x q identity matrix. 
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Let us assume that a point i is to be interpolated from q C neighbors. Define a 
composite grid consisting of the fine-grid in the interpolation region and progressively 
coarser grids as we move away from the target point i. We now use AMG on the 
composite grid with the current coarsening to compute an n x q matrix Q with the 
following property: 

Q minimizes p(Q T LQ) 

subject to 

(PQfPQ = I 

where p(A) denotes the spectral radius of A. Note that the q x q matrix PQ denotes 
the restriction of Q to the q interpolatory C points. 

The motivation here is to find the vectors that lie in the near null space of L 
that are very distinct on the C points. Included in the space determined by Q (i.e., 
its span) is the “eigenvector” p of L belonging to the smallest eigenvalue where the 
local normalization, ( Pp) T Pp = 1, is used. 

This scheme is reasonably cheap except that the usual coarse grids would 
introduce O(log n) complexity. McCormick notes that all of the coarse grids may not 
be needed for the interpolatory region and is hopeful that we’ll need just a few very 
coarse ones. In this manner, we can get the interpolation we want in a fairly general 
setting. 

If we simplify this local problem by interpreting the columns of Q to be vectors 
truly defined on the fine-grid then it should be fairly easy to implement. With only 
one target point, it’s easy to see. Then, Q is an n vector. If we now partition Q into 
its components u, belonging to the C points, and v belonging to all other points, then 
the constraint ( PQ) T PQ = I just becomes p T p = 1. 

If we partition L corresponding to the u-v partition, we have 

A B 
B T C 
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and it’s can be seen that the local problem is just the eigenvalue problem: 



minimize RQ(u) 



where 

u T Du 

T 

U 1 U 

and 

D = A-BC~ 1 B t . 

Note that this determines u, and then v can be computed by v = —C~ 1 B T u. Now, 
when we want more targets, we just solve for more eigenvectors of D. In fact, we 
want all q eigenvectors of the q x q matrix D. Once these eigenvectors are known, a 
least-squares procedure could be used to select interpolation weights, similar to the 
manner described in Section B. Preliminary experiments with the method outlined 
in this section have yielded some favorable indications, but there is much yet to be 
accomplished before the efficacy of the approach can be demonstrated. 

The three avenues of research outlined here are all active areas of effort. All 
are focused on the principle that accurate interpolation of “local” null-space and 
near null-space vectors will ensure accurate approximation of the global equivalents. 
Moreover, an increase in accuracy of interpolation for smooth errors results in a more 
robust and efficient V-cycle convergence making algebraic multigrid a competitive 
alternative to current methods. In fact, in some cases, there is no better matrix 
solver. Even though this research is in its infancy, indications now point toward 
eventual success in using these new ideas to produce a robust, accurate, and efficient 
interpolation. 
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GLOSSARY 



algebraic multigrid a numerical method used to solve particular algebraic (linear) 
systems with the “multigrid-style” cycling process 

boundary conditions a condition imposed on a bounding surface (in three dimen- 
sions) or line (in two dimensions) or at a bounding point (in one dimension) 
coarse-grid operator _A TO+1 = 7™ +1 A m /™ +1 

connections a point i £ Q, h is connected (directly) to a point j £ $l h if a 1 }- / 0 
diagonal dominance (strict) |a tt | > |aq| i = 1 , . . . ,n 

diagonal dominance (weak) |a;;| > \ a ij I i = 1 , . . . , n 

Dirichlet boundary conditions u = 0 on the boundary 
injection one type of a linear restriction operator 
iterative methods a.k.a relaxation 

Fourier modes vj = sin (j kir /N) where k is the wave number 
Galerkin condition A m+l = I™ +1 A m I™ + 15 known as the coarse-grid operator 
interpolation operator 7™ +1 : 7£ nm+1 — *■ H nm 
Laplace’s equation V 2 u = 0 

M-matrix a matrix which is positive definite and whose off-diagonal components 
are nonpositive 

Poisson’s equation V 2 u = 5, S ^ 0 
restriction operator 7™ +1 : TZ nm — > 7£ nm+1 

smoothing operator G : TZ n — > 7Z n , G acts on e to remove the oscillatory modes 
of the error 

smoothing property the property of an iterative method that quickly removes 
the oscillatory modes but leaves the smooth (or low-frequency) modes of the error 
strongly connected a point i is strongly connected to a point j if |a,j| ~ max|a t fc| 
for j £ C{ and k £ Ni 

variational properties 7™ +1 = (7™ +1 ) T and J 4 m+1 = I™ +1 A m I™ +1 
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wave number the kth mode consists of k/2 full sine waves on the interval (0, N) 
C- variables coarse level variables 
F- variables fine level variables 
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